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ABSTRACT 


During  the  three-year  period  of  this  project,  we  conducted  research  on  the  develop¬ 
ment,  analysis,  and  application  of  serial  and  parallel  adaptive  computational  strategies  for 
solving  transient  and  steady  partial  differential  systems.  Concentrating  on  high-order 
methods  and  adaptive  approaches  that  combine  mesh  refinement  and  coarsening  (h- 
refinement),  order  variation  (p-refinement).  and,  occasionally,  mesh  motion  (r-refineraent), 
we  addressed  problems  in  combustion,  materials  science,  and  compressible  fluid  mechan¬ 
ics.  Special  spatially-discrete  finite  element  Galerkin  methods  were  considered  for  the 
parallel  and  adaptive  solution  of  hyperbolic  conservation  laws.  Improved  solution-limiting 
and  error-estimation  strategies  increased  the  accuracy  and  efficiency  of  these  methods 
which  are  being  applied  to  two-  and  three-dimensional  compressible  flow  problems. 
Adaptive  techniques  for  dissipative  systems  are  being  applied  to  problems  m  the  manufac¬ 
ture  of  ceramic  composite  media. 
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1.  Research  Status  and  Key  Results. 

Research  on  this  project  involved  the  development,  analysis,  and  application  of  adap¬ 
tive  finite  element  methods  for  the  numerical  solution  of  transient  and  steady  partial 
differential  systems.  As  a  point  of  reference,  consider  th  vector  partial  differential  system 

Mu,  +  f(x,/,u,Vu)  =  Vd(x,r,u,Vu),  xe  Q,  t  >  0.  (1) 

where  t  denotes  time,  x  denotes  spatial  position,  V  is  the  spatial  gradient  or  divergence 
operator,  Q  is  a  bounded  one-,  two-,  or  three-dimensional  region,  M  is  a  positive  semi- 
definite  ma.ss  matrix,  and  the  functions  f  and  d  describe  the  effects  of  reaction,  convection, 
and  diffusion.  Specific  techniques  aim  to  solve  parabolic,  elliptic,  and  hyperbolic  prob¬ 
lems  that  are  obtained  from  ( 1 )  by  suitable  restrictions  of  M  and  d. 

Our  concentration  involved  high-order  methods  and  adaptive  approaches  that  com¬ 
bined  mesh  refinement  or  coarsening  (h-refinement),  order  variation  (p-refinement),  and, 
occasionally,  mesh  motion  (r-refinement).  The  particular  combination  of  h-  and  p- 
refinement  typically  produces  (hp-refinement)  schemes  with  exponential  convergence  rates. 
Temporal  integration  relied  on  either  a  method  of  lines  (MOL)  framework  where  local 
spatial  enrichment  is  combined  with  global  temporal  enrichment  or  a  local  refinement 
method  (LRM)  which  features  local  spatial  and  temporal  enrichment.  Our  work  on  paral¬ 
lel  adaptive  strategies  also  concentrated  on  high-order  techniques  with  the  goal  of  develop¬ 
ing  parallel  hp-refinement  strategies.  Given  their  superior  performance  on  serial  comput¬ 
ers,  consideration  of  other  sub-optimal  strategies  would  be  dubious. 

1.1.  Method  of  Lines  Techniques. 

Adaptive  MOL  procedures  utilize  the  power  of  sophisticated  ordinary  differential 
equations  software  and  avoid  the  need  to  develop  special  techniques  for  time  integration. 
Using  the  differential-algebraic  systems  software  DASSL  developed  by  Petzold  and  the 
DASSL  implementation  within  the  SPRINT  package,*  we  studied  adaptive  finite  element 
strategies  and  developed  software  for  parabolic  systems  that  combined  h-,  p-,  and,  in  one 
dimension,  r-refinement  [8,  11,  12,  16,  191.^  The  two-dimensional  software  uses  a  finite 
quadtree  structure  to  manage  both  mesh  generation  and  refinement  An  initial  unstructured 
mesh  of  triangular  and/or  quadrilateral  elements,  created  from  terminal  quadrants  of  the 
tree  structure,  may  be  further  refined  or  coarsened  by  adding  or  deleting  quadrants  to  or 
from  the  tree  using  solution-based  h-refinement  indicators.  By  modifying  quadrants 
instead  of  the  mesh,  we  maintain  an  underlying  uniform  structure  while  performing 

*  M.  Berzins  and  RM.  Furzeiand,  ‘‘A  User’s  Manual  for  SPRINT  -  A  Versatile  Software  Package  for 
Solving  Systems  of  Algebraic,  Ordinary  and  Partial  Differential  Equations:  Part  1  -  Algebraic  and  Ordinary 
Differential  Equations,”  Shell  Research  Ltd.,  Thornton  Research  Centre,  Amsterdam.  1985. 

^  Citations  correspond  to  the  publications  and  manuscripts  listed  in  Section  3. 
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computations  on  an  unstructured  mesh. 

Estimates  of  spatial  discretization  errors  are  obtained  by  a  p-retinement  technique 
using  piecewise  polynomial  approximations  of  one  degree  higher  than  those  used  for  the 
solution  to  solve  local  problems  with  forcing  proportional  to  elemental  or  edge  residuals. 
Errors  of  odd-degree  finite  element  solutions  are  approximated  using  edge  residuals  while 
neglecting  elemental  residuals.  Those  of  even  degree  solutions  follow  the  opposite  course. 
Error  estimates  computed  in  this  manner  converge  in  energy  at  the  correct  rate  to  the 
actual  finite  element  spatial  errors  under  h-refinement  [11].  Our  results  show  that  this 
edge  or  elemental  “superconvergence”  is  robust  in  time;  hence,  temporal  variation  of  spa¬ 
tial  errors  may  be  neglected  [11,  12].  This  enables  us  to  determine  error  estimates  by 
solving  local  elliptic  instead  of  local  parabolic  problems.  The  advantage  of  this  approach 

is  that  error  estimates  may  be  computed  only  when  needed  and  not  at  every  time  step 

selected  by  the  ordinary  differential  equations  software.  In  typical  applications  [8,  11,  12], 
this  reduces  the  time  required  to  obtain  an  error  estimate  by  approximately  thirty  percent. 

1.1.1.  Case  Study:  Boundary  Resonance. 

The  problem 

Ut  +  xux  =  tu„,  -1  <  jr  <  1,  r  >  0,  (2a) 

«(-!,/)  =  3,  «(l.r)  =  5.  r>0,  (2b) 

u{x,0)  =  x+4,  -l<x<l.  (2c) 

is  an  example  of  boundary  resonance^  in  the  singularly-perturbed  limit  0  <  e  ^  1.  As  t 
increases,  the  smooth  initial  data  evolves  to  the  steady  state 

u(x)  =  xe^^'^M  (1,-^,--^)  +  4  (2d) 

2  2£ 

where  is  Kummer’s  function.  This  steady  solution  has  boundary  layers  of 

width  0(Ve)  at  both  endpoints.  The  problem  (2),  connected  with  e.xit-time  probability 
determination,  is  challenging  because  (i)  singular-perturbation  theory  does  not  suffice  to 
determine  the  interior  solution  (i.e.,  the  solution  on  -1  <x  <  1)  and  (ii)  the  transient  solu¬ 
tion  evolves  to  a  steady  state  at  an  exponentially  slow  rate  (i.e.,  in  time). 

Adjerid  et  al.  [12,  19]  computed  solutions  of  (2)  using  adaptive  h-,  hr-,  and  hp- 
refinement  with  e  =  0.001  and  several  error  tolerances  in  H^.  Integration,  starting  with  a 


^  C.  Holland,  Personal  communication,  1989. 
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ten-eletnent  mesh,  was  continued  to  t  =  20,  which  was  taken  as  steady  state.  The  solution 
generated  by  adaptive  hp-refinement  with  an  error  tolerance  of  0.03  relative  to  ||u||i  is 
shown  at  r  =  0  and  20  in  Figure  1. 


1 

55 


Figure  1.  Finite  element  solution  of  (2)  at  r  =  0  and  20  using  adaptive  hp- 
refinement  with  a  relative  tolerance  of  0.03  on  a  10-element  initial  mesh  (left). 
Relative  error  as  a  function  of  work  (O  using  adaptive  h-,  hr-,  and  hp-refinement 
(right). 


The  discretization  error  is  shown  as  a  function  of  the  total  work  [12]  Cf)  for  solutions 
computed  by  adaptive  h-,  hr-,  and  hp-refinement  at  the  right  of  Figure  1.  Accuracy  of  all 
methods  are  comparable  for  the  larger  errors,  but  hp-refinement  is  much  more  efficient  at 
higher  accuracy.  A  “warm  restarting”  procedure  [8]  enables  the  backward  difference 
software  within  DASSL  to  resume  the  temporal  integration  at  a  high  order  whenever  h- 
refinement  has  been  performed  and  this  has  improved  the  performance  of  the  b-  and  hp- 
refinement  schemes  and  has  reduced  the  need  for  mesh  motion. 

Consider  a  related  two-dimensional  problem'* 


^  C.  Holland,  Personal  communication,  1991. 
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u,  ■¥({x,yyVu  =zAu,  (x,y)eQ,  t>0,  (3a) 

with 

f(jc,v)  =  [-V  -  .r(l  -  r  \  X  -  v(l  -  r^)]^.  (3b) 

Q  =  { (jr,y)  I  -2  <  <  i.5,  -3  <  v  <  3  },  and  =  x^  +  y^.  With  this  resonance  problem, 

the  solution  in  the  interior  of  the  domain  is  determined  by  the  data  at  the  point(s)  on  the 
boundary  that  is  (are)  closest  to  the  origin. 

The  solution  shown  in  Figure  2  was  obtained  by  integrating  the  transient  parabolic 
problem  with  e  =  0.01  to  steady  state  using  hp-refinement  with  a  relative  error  tolerance  of 
0.02  in  Mesh  refinement  with  relatively  low-order  methods  occurs  in  boundary  layers 
near  y  =  ±3  and  x  =2.  The  smooth  solution  in  the  interior  and  near  or  =  1.5  was 
obtained  on  a  coarser  mesh  principally  with  p-refinement.  Standard  numerical  techniques 
for  solving  this  problem  would  produce  non-physical  oscillations  unless  the  mesh  was 
sufficiently  fine  within  the  boundary  layers  but,  with  hp-refinement,  the  solution  has  been 
captured  correctly  in  near  optimal  fashion. 

As  a  final  example  of  resonance,  consider  the  elliptic  problem^ 

xUx+yUy-eAu,  (Ar,y)eQ.  (4a) 

with  Q  =  { (x  ,y )  1  -2  <  X  <  I,  -3  <  y  <  3  Problems  (2)  and  (3)  were  solved  using 
adaptive  techniques  with  standard  finite  element-Galerlcin  methods.  In  order  to  reduce  the 
dependence  on  proper  mesh  placement  or  sophisticated  temporal  embedding,  we  have 
developed  special  quadrature  rules  that  produce  stable  numerical  strategies  for  a  variety  of 
singularly-perturbed  problems.  This  technique  will  be  described  in  a  forthcoming 
manuscript;^  however,  for  the  moment,  we  note  that  these  rules  are  approximately  Radau 
quadrature  for  this  example.  The  solution  of  (4)  shown  in  Figure  3  was  obtained  with 
e  =  0.001  using  adaptive  h-refinement  on  a  20x  20  initial  mesh  with  piecewise-linear 
approximations  and  the  special  quadrature  rules.  Boundary  layers  on  the  three  edges  that 
are  farthest  from  the  origin  are  sharp,  no  spurious  oscillations  are  present,  and  the  correct 
interior  solution  is  obtained. 

These  resonance  problems  are  extremely  difficult  Convergence  to  steady  state  at  an 
exponentially  slow  rate  implies  that  discrete  versions  of  steady  problems,  such  as  (4),  are 
extremely  ill  conditioned.  Standard  and  adaptive  techniques  must  fail  when  e  is 
sufficiently  small.  The  adaptive  h-  and  hp-refinement  strategies  have  enabled  us  to 

^  C.  Holland,  Personal  conununication,  1989. 

^  S.  Adjerid,  M.  Aiffa,  and  J.E.  Flaherty.  “Adaptive  Numerical  Procedures  for  Singularly-Perturbed  Par¬ 
tial  Differential  Equations,"  Workshop  on  Perturbation  Methods  in  Physical  Mathematics,  Rensselaer  Po¬ 
lytechnic  Institute,  Troy,  June  23-26,  1993.  Manuscript  in  preparation. 


Figure  2.  Mesh  (left)  and  steady  solution  (right)  of  (3)  with  e  =  0.01  using 
adaptive  hp-refinement. 

produce  solutions  at  smaller  values  of  £  than  would  otherwise  be  possible.  The  special 
quadrature  rules  eliminate  all  oscillations.  However,  additional  problem-dependent  infor¬ 
mation^  would  be  necessary  in  order  to  obtain  solutions  for  very  small  values  of  e. 

1.1.2.  Case  Study:  Materials  Processing. 

Working  with  Rensselaer  Materials  Scientists  William  Hillig,  John  Hudson,  and  Nag 
Patibandla,  we  have  been  using  our  adaptive  software  to  study  fabrication  problems 

*  cf.,  e.g.,  M.  Ward,  "On  Exponentially  Slow  Internal  Layer  Motion,"  Workshop  on  Perturbation 
Methods  in  Physical  Mathematics,  Rensselaer  Polytechnic  Institute,  Troy,  June  23-26,  1993. 
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Figure  3.  Mesh  (left)  and  solution  (right)  of  (4)  with  e  =  0.001  using  adaptive 
h-refinement  with  piecewise-linear  approximations  on  a  20  x  20  initial  mesh  and 
special  quadrature  rules. 


associated  with  ceramic  composites.  The  goal  is  to  manufacture  high-temperature  compo¬ 
sites  having  alumina  (-4/20  3)  fibers  and  a  molybdenum  di-silicide  {MoSi^  matrix.  As  a 
preliminary  study,  we  have  been  modeling  the  molybdenum  di-silicide  fabrication  by  a 
“reactive  vapor  infiltration”  (RVD  process  developed  by  Hillig  and  Patibandla.  Using 
this  process,  molybdenum  powder  is  pressed  (to  45%  porosity)  into  a  pellet  which  is 
exposed  at  approximately  1200®  C  to  a  5/C/4  and  H2  flow.  Silicon  gas  is  produced  on 
the  pellet’s  surface.  It  infiltrates  the  pellet  and  reacts  with  the  molybdenum  to  form  sili- 
cide  layers  of  MosSi^  and  M0S12.  Pores  between  particles  are  filled  by  volume  expansion 
as  the  silicides  are  produced  and,  as  an  initial  approximation,  we  assumed  that  solid-state 
diffusion  is  the  dominant  means  of  infiltration.  With  this  assumption,  considerations  of 
mass  conservation  imply 

</©  + j^J,  ds  =  /J/p,  </o),  I  =1,  2, —.W, 


(5a) 


(5b) 
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N  N 

=  1,  2p,  =0. 

(=1  isl 

where  V  ,  J, .  and  p,  are,  respectively,  the  mass  fraction,  flux,  and  production  rate  of 
spcL  wi  /  ,  /■  =  1,  2,  •••,  N,  and  p  is  the  mixture  density.  Assuming  Fickian  diffusion  and 
first-order  reactions  yields  the  dimensionless  system  for  the  silicon-molybdenum  reaction 


3m  I  dM 5  Ml 

*  V7  r\  r7  -  v  ^  n  r7  ^  ^ 


dt 

-  VD,Vm,  = 

-(3w'  1  +  7kh’2), 

dt 

-  V  D2VU2  =  5^^  IOV2’ 

(6a,b) 

dM3 

-  = 

M3 

3m  4 

M4 

dt 

(h',-k:w2), 

M, 

dt 

-  V  D^Vu4  =  -5  w, 

M , 

(6c,d) 

4 

5^M,=p,  H’,=M,M4,  H’2  =  M,M3.  (6e,f.g) 

1=1 

The  volume  fraction  u, ,  diffusivity  D, ,  and  molar  mass  M,  of  species  / ,  /  =1,2,3.  4, 
respectively,  correspond  to  those  of  Si,  MoSi2,  Mo^Si^,  and  Mo.  The  variable  k 
identifies  the  ratio  of  the  production  rate  of  Mo5i2  to  that  of  Mo^Si^  P  is  the  dimen¬ 
sionless  mixture  density. 

We  used  this  model  with  our  adaptive  software  to  solve  several  one-  and  two- 
dimensional  problems.  Results  using  adaptive  h-refinement  shown  in  Figures  4  and  5, 
respectively,  exhibit  the  thickness  of  the  MoSi2  layer  at  three  times  and  the  concentrations 
of  Si,  MoSii,  Mo^Si^,  and  Mo  after  3.5  hours  of  reaction.  In  each  drawing,  high  concen¬ 
trations  are  colored  red  and  low  ones  arc  colored  blue.  The  mesh  used  for  the  computa¬ 
tion  at  the  indicated  time  has  been  superimposed.  The  bottom  and  right  edges  of  each 
figure  are  lines  of  symmetry.  The  silicide  layers  progress  inward  from  the  top  and  left 
edges  to  build  the  MOSi2  Mo^Si^  layers.  Reaction  fronts  are  sharp  and  the  mesh  is 

concentrated  in  and  following  these  regions.  The  Mo^Si^  layer  is  a  narrow  zone  that  pre¬ 
cedes  the  MoSi2  layer  into  the  unreacted  molybdenum. 

The  model  and  its  adaptive  solution  correctly  predicted  the  thickness  and  growth  rate 
of  the  MoSi2  layer  to  approximately  5%  of  measured  values.  Volume  expansion  also 
agreed  with  experimental  observations.  The  generality  of  our  software  makes  it  easy  to 
change  models  and  include  more  complex  and  realistic  effects.  Recently,  we  were  able  to 
show  that  the  production  times  of  MoSi2  could  be  reduced  (potentially  halved)  by  starting 
with  a  powder  composed  of  a  mixture  of  Mo,  MoSi2,  and/or  Mo^Siy  Hillig’s  and 
Patibandla’s  initial  experiments  are  confirming  these  computational  predictions.  Since  pro¬ 
duction  times  are  approximately  one  day  for  realistic  layer  thicknesses,  a  factor  of  two 
reduction  is  very  significant. 
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1.2.  Local  Refinement  Methofls. 

1.2.1.  Parabolic  Systems. 

While  MOL  procedures  greatly  simplify  temporal  integration,  there  are  many  situa¬ 
tions  where  they  would  have  sub-optimal  performance  and,  for  this  reason,  we  have  also 
been  studying  LRMs.  A  posteriori  error  estimates  may  be  used  to  schedule  different  tem¬ 
poral  and  spatial  enrichment  strategies  in  different  regions  of  the  spatial  domain  with  the 
possibility  of  providing  large  efficiency  improvements.  Given  the  need  to  interrupt  the 
temporal  integration  for  h-refinement  more  often  than  with  a  MOL  procedure,  one-step 
integration  methods  may  be  preferable  to  the  backward  difference  methods  that  are  com¬ 
monly  used  with  MOLs  [2,  4|.  With  this  possibility,  we  have  developed  a  variant  of  the 
singly  implicit  Runge-Kutta  (SIRK)  method  that  offers  several  advantages  when  used  in 
conjunction  with  hp-refinement  methods  for  parabolic  problems  [18].  SIRKs  have  high- 
stage  order,  A-stability  for  orders  one  through  eight,  efficiencies  close  to  those  of  back¬ 
ward  difference  software,  embedded  error  estimates,  and  a  locality  that  is  suitable  for 
parallel  computation.  Temporal  error  estimates  of  each  SIRK  stage  may  be  computed  for 
the  cost  of  one  additional  function  evaluation  (18).  This  permits  the  acceptance  of  solu¬ 
tions  at  partial  time  steps  whenever  accuracy  is  not  satisfactory  at  the  final  time  [25],  An 
adaptive  hp-refinement  strategy  based  on  SIRKs  and  hierarchical  finite  element  approxima¬ 
tions  is,  to  our  knowledge,  the  first  instance  where  hp-refinement  has  been  used  in  both 
space  and  time  [25],  Several  enrichment  strategies  were  developed  and  the  most  success¬ 
ful  of  these  are  being  used  to  develop  a  three-dimensional  solution  package. 

1.2.2.  Hyperbolic  Systems. 

While  hp-refinement  is  reasonably  well  understood  for  elliptic  problems,  this  is  far 
from  true  for  time-dependent  problems,  particularly  hyperbolic  systems.  Extensions  of 
finite  difference  methods  to  higher  order  involve  enlarging  their  computational  stencils, 
which  would  be  unsuitable  for  both  parallel  and  unstructured-mesh  computation.  Cock- 
bum  and  Shu  (cf.,  e.g.,  [14,  15,  22,  24])  combine  a  local  basis  of  Legendre  polynomials 
with  a  solution  limiting  scheme  to  create  a  method  that  has  both  high  order  and  a  compact 
stencil.  In  his  Rensselaer  dissertation,  Biswas^  used  this  scheme  to  develop  adaptive  h- 
refinement  methods  for  one-  and  two-dimensional  problems.  Subsequent  extensions 
include  improved  projection  limiting  schemes  near  discontinuities,  spatial  a  posteriori  error 
estimation  schemes  based  on  p-refinement,  adaptive  p-refinement,  and  parallel  solution 
sffategies  [15,  22,  24],  Spatial  error  estimates  were  computed  using  Radau  quadrature 
points  as  described  in  Section  1.1.1. 


^  Cf.  the  list  of  dissertations  in  Section  4. 
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1.3.  Parallel  Computation. 

The  Construction  of  parallel  adaptive  techniques  is  a  challenge  since  parallelism  and 
adaptivity  are  at  odds.  Adaptivity  employs  complex  heuristic  decision  processes  that 
introduce  serial  bottlenecks  and.  hence,  complicate  parallelization.  This  situation  is  further 
exacerbated  when  unstructured  grids  are  used  as  part  of  the  computation.  Nevertheless,  as 
noted,  parallelization  of  these  most  efficient  techniques  must  be  explored  since  paralleliz¬ 
ing  sub-optimal  strategies  will  result  in  lower  performance. 

1.3.1.  Elliptic  Systems. 

The  need  for  synchronization  on  shared-memory  computers  can  be  reduced  by  con¬ 
structing  a  list  of  elements  that  have  independent  basis  support  that  may  be  assembled  and 
solved  in  parallel.  The  task  of  separating  or  “coloring”  the  elements  should  be  (i)  com¬ 
putationally  efficient  and  (ii)  employ  a  small  number  of  colors  to  maintain  granularity. 
Algorithms  to  color  arbitrary  finite  element  meshes  with  a  minimum  number  of  colors  are 
data  and  manipulation  intensive;  however,  meshes  produced  by  finite  quadtree  and  finite 
octree  procedures  have  a  regular  quadrant  or  octant  structure  that  may  be  used  to  advan¬ 
tage.  Thus,  we  color  the  quadrants  or  octants  of  the  tree  structure  instead  of  the  finite  ele¬ 
ments  and  thereby  simplify  the  coloring  process  at  the  expense  of  creating  a  minor  load 
imbalance  due  to  the  variation  in  the  number  of  elements  per  nuadrant  or  octant. 

Benantar  et  al.  [3,  6,  7.  9,  10,  17]  described  a  six-color  procedure  for  coloring  quad¬ 
tree  structures  that  was  used  to  solve  elliptic  problems  with  a  piecewise-linear  basis.  Con 
jugate  gradient  iteration  with  element-by-element  and  symmetric  successive  over  relaxation 
(SSOR)  preconditionings  were  used  to  solve  the  resulting  linear  algebraic  systems.  The 
procedures  worked  well  with  the  SSOR  prcconditioner  providing  the  best  performance; 
however,  performance  degraded  when  high-order  approximations  and  adaptive  p- 
refinement  were  added  to  the  solution  process.  Since  hierarchical  bases  of  degree  greater 
than  unity  involve  edge-  and  element-level  interactions,  an  edge  coloring  procedure  should 
improve  parallel  performance. 

Benantar,  in  his  Rensselaer  dissertation,  considered  edge  coloring  procedures  for 
two-dimensional  meshes  of  triangular  elements.  Creating  a  new  discrete  structure,  called  a 
triangle  graph,  he  showed  that  mesh  edges  could  be  colored  using  a  maximum  of  three 
colors  by  a  procedure  that  had  at  most  quadratic  time  complexity  [17,  21],  Contrary  to 
the  quadrant-coloring  procedure,  performance  improves  with  increasing  polynomial  degree; 
thus,  this  is  an  effective  procedure  for  parallel  p-refinemenL 

1.3.2.  Hyperbolic  Systems. 

We  are  using  the  spatially  discontinuous  Galerkin- finite  element  procedure  described 
in  Section  1 .2.2  to  construct  parallel  h-  and  p-refinement  techniques  on  MIMD  systems  for 
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two-  and  three-dimensional  conservation  laws.  The  two-dimensional  p-refinement 
software  uses  a  tiling  strategy  developed  by  Wheat  [24]  to  appraise  neighborhood  work¬ 
loads  and  migrate  elements  to  nearby  processors  as  appropriate.  Migration  is  far  more 
efficient  with  adaptive  computation  than  global  optimization  methods.  Computational 
experiments  using  a  ten-dimensional  NCube  2  hypercube  at  Sandia  National  Laboratories 
indicated  a  52  percent  increase  in  performance  with  this  tiling  algorithm  relative  to  no 
migration.  We  expect  further  performance  improvements  to  occur  by  predicting  future 
workload  changes  using  solution  and  error  information. 

Extending  the  spatially  discontinuous  Galerkin  procedure  to  three  dimensions,  we 
consider  parallel  adaptive  h-refinement  on  octree-structured  grids.  An  example  involving 
the  supersonic  (Mach  two)  flow  over  a  cone  is  shown  in  Figures  6  and  7.  In  Figure  6,  we 
show  the  initial  and  final  meshes  used  for  this  computation.  Refined  meshes  were  created 
by  introducing  additional  octants  into  the  tree  structure  and  generating  a  new  mesh  locally. 
Steady  solutions  were  obtained  by  integrating  the  transient  Euler  equations  to  steady  state 
with  local  time  steps  proportional  to  tree  level.  Load  balancing  on  a  CM-5  Connection 
Machine  was  performed  by  a  traversal  of  the  octree  structure  using  a  work  metric  at  each 
tree  node.  Pressure  contours  and  the  shock  surface,  shown  in  Figure  7,  indicate  that  the 
mesh  is  being  refined  in  the  vicinity  of  the  leading  bow  shock.  Higher-order  methods  and 
p-refinement  will  be  added  to  the  method  shortly. 
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Figure  6.  Coarse  (top)  and  fine  (bottom)  meshes  used  to  compute  the  Mach 
two  Row  past  -  cone  having  a  half  angle  of  10'^. 
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Figure  7.  Shock  surface  and  pressure  contours  found  when  computing  the 
Mach  two  flow  past  a  cone  having  a  half  angle  of  10'^. 


-  17  - 


2.  Interactions. 

Key  invited  lectures  by  J.  E.  Flaherty  on  material  connected  with  the  research  sup¬ 
ported  by  this  grant  are  listed  below.  (Invited  university  colloquia  are  only  presented  for 

the  last  twelve  months.) 

1.  “Adaptive  Methods  and  Parallel  Computation  for  Panial  Differential  Equations,” 
Special  Session  on  Adaptive  Methods  at  the  Eighth  Army  Conference  on  Applied 
Mathematics  and  Computing,  Cornell  University,  Ithaca,  June  19-22.  1990. 

2.  “Symbolic  and  Parallel  Adaptive  Methods  for  Partial  Differential  Equations,” 
Workshop  on  the  Integration  of  Numerical  and  Symbolic  Computing  Methods,  Sara¬ 
toga  Springs,  July  8-11,  1990. 

3.  “Advances  in  Parallel  Processing  for  Electromagnetic  Field  Applications,”  Fourth 
Biennial  IEEE  Conference  on  Electromagnetic  Field  Computation,  Toronto.  October 
22-24,  1990. 

4.  “Adaptive  Finite  Element  Schemes  for  Parabolic  Partial  Differential  Systems  with 
H-,  P-,  and  R-Refinement,”  AJCHE  Annual  Meeting,  Chicago,  November  12-16, 
1990. 

5.  “Combined  Symbolic  and  Numerical  Adaptive  Methods  for  Partial  Differential  Equa¬ 
tions,”  Workshop  on  Problem  Solving  Environments,  National  Science  Foundation, 
Washington,  April  11-12,  1991. 

6.  “Adaptive  Methods  for  PDEs,”  Conference  on  Experimental  Mathematics:  Computa¬ 
tional  Issues  in  Nonlinear  Science,  Los  Alamos  National  Laboratory,  Los  Alamos, 
May  20-24,  1991. 

7.  “Adaptive  h-  and  p-Refinement  Finite  Element  Methods  for  Singularly  Perturbed 
Parabolic  Systems,”  ICIAM  91,  Minisymposium  on  "Numerical  Methods  for 
Singularly-Perturbed  Problems,"  Washington,  July  8-12,  1991. 

8.  “A  Posteriori  Error  Estimation  for  Parabolic  Systems”  and  “High-Order  Adaptive 
Methods  for  Parabolic  Systems,”  First  U.  S.  National  Congress  on  Computational 
Mechanics,  Chicago,  July  21-24,  1991. 

9.  “Parallel  Adaptive  Methods  for  Partial  Differential  Equations,”  Workshop  on  Paral¬ 
lel  Scientific  Computation,  Rensselaer  Polytechnic  Institute,  Troy  October  31, 
November  1,  1991. 

10.  “Adaptive  Method-of-Lines  Techniques  for  Parabolic  Systems,”  First  International 
Symposium  on  Method  of  Lines,  Surfaces  and  Dimensional  Reduction  in  Computa¬ 
tional  Mathematics  and  Mechanics,  Athens,  November  12-16,  1991, 
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11.  “High-Order  Adaptive  Methods  for  parabolic  and  Elliptic  Systems,”  New  Jersey 
Institute  of  Technology,  Newark,  April  20,  1992. 

12.  “Parallel  Adaptive  Techniques,”  Workshop  on  Adaptive  Methods  for  Partial 
Differential  Equations,  Rensselaer  Polytechnic  Institute,  Troy,  May  17-20,  1992. 

13.  “Parallel  High-Order  Finite  Element  Methods  for  Conservation  Laws,”  /MACS 
Seventh  Internation  Conference  on  Computer  Methods  for  Partial  Differential  Equa¬ 
tion  "i,”  New  Brunswick,  June  22-24,  1992. 

14.  “Adaptive  Methods  for  Time- Dependent  Partial  Differential  Equations,”  SIAM  40th 
Anniversary  Meeting,  Los  Angeles,  July  20-24  1992. 

15.  “Adaptive  Methods  for  Time-Dependent  Partial  Differential  Equations,”  ICASE, 
NASA  Langley  Research  Center,  Hampton,  July  23,  1992. 

16.  “Adaptive  Methods  for  T'me-Dependent  partial  Differential  Equations.  Part  I:  Basic 
Strategies  and  Error  Estimation”  and  “Part  II;  High-Order  Methods  and  Parallel 
Computation,”  Dutch  Numerical  Analysis  Seminar,  Woudschoten  Conference  Center. 
Zeist,  October  5-7,  1992. 

17.  “High-Order  Adaptive  Methods  for  Time-Dependent  Partial  Differential  Equations,” 
Washington  University,  St.  Louis,  October  29,  1992. 

18.  “Coloring  Procedures  for  Finite  Element  Computation  on  Shared-Memory  Comput¬ 
ers,”  Symposium  on  Adaptive,  Multilevel,  and  Hierarchical  Computational  Stra¬ 
tegies,  ASME  Winter  Annual  Meeting,  Anaheim,  November  8-13,  1992. 

19  “Adaptive  and  Parallel  High-Order  Methods  for  Conservation  Laws,”  Mississippi 
State  University,  Mississippi  State,  January  22,  1993. 

20.  “Adaptive  and  Parallel  High-Order  Methods  for  Conservation  Laws,”  University  of 
Texas,  Austin,  March  13-22,  1993. 

21.  “Adaptive  and  Parallel  High-Order  Methods  for  Conservation  Laws,”  University  of 
Kentucky,  Lexington,  April  21,  1993. 

22.  “Adaptive  hp-Refinement  Methods  for  Parabolic  Partial  Differential  Equations,” 
1993  AFOSR  Workshop  on  Computational  Mathematics,  Washington  University,  St 
Louis,  May  19-21,  1993. 

Other  interactions  follow. 

1.  J.  E.  Flaherty  and  B.  K.  Szymanski  organized  a  workshop  on  Parallel  Scientific  Com¬ 
putation  which  was  held  at  the  Rensselaer  Polytechnic  Institute,  October  31  and 
November  1,  1991. 
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2.  J.  E.  Flaherty  and  M.  S.  Shephard,  organized  the  Workshop  on  Adaptive  Methods  for 
Partial  Differential  Equations  which  was  held  at  the  Rensselaer  Polytechnic  Institute, 
Troy,  May  17-20,  1992. 

3.  J.  E.  Flaherty  organized  a  special  session  on  Adaptive  Methods  for  Time-Dependent 
and  Nonlinear  Partial  Differential  Systems  which  was  held  at  the  I  MACS  Seventh 
Intemation  Conference  on  Computer  Methods  for  Partial  Differential  Equations, 
New  Brunswick,  June  22-24,  1992. 

4.  Lenart  Johnson  of  Thinking  Machines  Corporation  visited  J.  E.  Flaherty  at  Rensselaer 
on  September  10,  1992.  Flaherty  and  graduate  student  Raymond  Loy  visited  John¬ 
son,  Zdenek  Johann,  et  al.  at  Thinking  Machines  on  March  3,  1993.  Both  visits 
involved  discussions  of  parallel  adaptive  techniques  for  partial  differential  equations. 

5.  L.  R.  Petzold  of  the  University  of  Minnesota  visited  J.  E.  Flaherty  on  September  30, 
1992  to  present  a  colloquium  and  discuss  research  on  parallel  adaptive  techniques  for 
partial  differential  equations. 

6.  J.  E.  Flaherty  attended  the  AFOSR  Electronic  Prototyping  Initiative  Review  Meeting 
at  the  University  of  Michigan,  Ann  Arbor,  May  27,  1993. 
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3.  Publications  and  Manuscripts 

Publications 

1.  PC.  Amey  and  J.E.  Flaherty,  “An  Adaptive  Mesh  Moving  and  Local  Refinement 
Method  for  Time-Dependent  Partial  Differential  Equations,"  ACM  Trans.  Math. 
Software,  16  (1990),  pp.  48-71. 

2.  P.K.  Moore  and  J.E.  Flaherty,  “A  Local  Refinement  Finite  Element  Method  for 
One-Dimensional  Parabolic  Systems,”  SIAM  J.  Numer.  Anal,  27  (1990),  pp.  1422- 
1444. 

3.  M.  Benantar,  R.  Biswas,  J.E.  Raherty,  and  M.  S.  Shephard,  “Parallel  Computation 
with  Adaptive  Methods  for  Elliptic  and  Hyperbolic  Systems,”  Comp.  Meths.  Appl. 
Mech.  and  Engrg.,  82  (1990),  pp.  73-93. 

4.  P.K.  Moore  and  J.E.  Raherty,  “Adaptive  Local  Overlapping  Grid  Methods  for  Para¬ 
bolic  Systems  in  Two  Space  Dimensions,”  J.  Comp.  Phys.,  98  (1992),  pp.  54-63. 

5.  D.C.  Amey,  R.  Biswas,  and  J.E.  Raherty,  “An  Adaptive  Mesh  Moving  and 
Refinement  Procedure  for  One-Dimensional  Conservation  Laws,”  Appl  Numer. 
Math.,  11  (1993),  pp.  259-282. 

6.  M.  Benantar  and  J.E.  Raherty,  “A  Six-Color  Procedure  for  the  Parallel  Solution  of 
Elliptic  Systems  Using  the  Finite  Quadtuw;  Structures,”  in  J.  Dongarra,  P.  Messing, 
D.C.  Sorensen,  and  R.G.  Voigt,  Eds.,  Proc.  Fourth  SIAM  Conf.  Parallel  Proc.  for 
Sci.  Comput.,  23  (1990),  pp.  231-236. 

7.  J.E.  Raherty,  M.  Benantar,  R.  Biswas,  and  P.K.  Moore,  “Symbolic  and  Parallel 
Adaptive  Methods  for  Partial  Differential  Equations,”  Chap.  12  in  B.R.  Donald,  D. 
Kapur,  and  J.L.  Mundy,  Eds.,  Symbolic  and  Numerical  Computation  for  Artificial 
Intelligence,  Academic  Press,  London,  pp.  277-302,  1992. 

8.  J.E.  Raherty  and  Y.J.  Wang,  “Experiments  with  an  Adaptive  h-,  p-,  and  r- 
Refinement  Finite  Element  Method  for  Parabolic  Systems,”  in  G.D.  Byme  and  W.E. 
Schiesser,  Eds.,  Recent  Developments  in  Numerical  Methods  and  Software  for 
ODEslDAEslPDEs,  World  Scientific,  Singapore,  pp.  55-80,  1992. 

9.  R.  Biswas,  M.  Benantar,  and  J.E.  Raherty,  “Adaptive  Methods  and  Parallel  Compu¬ 
tation  for  Partial  Differential  Equations,”  Trans.  Eighth  Army  Conf.  on  Applied 
Maths,  and  Comput.,  ARO  Report  91-1,  U.S.  Array  Research  Office,  Research  Trian¬ 
gle  Park,  pp.  531-552,  1991. 

10.  M.  Benantar,  R.  Biswas,  and  J.E.  Raherty,  “Advances  in  Adaptive  Parallel  Process¬ 
ing  for  Field  Applications,”  IEEE  Trans.  Magnetics,  27  (1991),  pp.  3768-3773. 
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11.  S.  Adjerid,  J.E.  Raherty,  and  Y.J.  Wang,  “A  Posteriori  Error  Estimation  with  Finite 
Element  Methods  of  Lines  for  One-Dimensional  Parabolic  Systems,”  Numer.  Math., 
65  (1993),  pp.  1-21. 

12.  S.  Adjerid,  J.E.  Raherty,  P.K.  Moore,  and  Y.J.  Wang,  ‘‘High-Order  Adaptive 
Methods  for  Parabolic  Systems,”  Physica  D.  60  (1992),  pp.  94-11 1. 

13.  B.E.  Webster,  M.S.  Shephard,  Z.  Rusak,  and  J.E.  Raherty,  ‘‘An  Adaptive  Finite  Ele¬ 
ment  Methodology  for  Rotorcraft  Aerodynamics,”  in  R.  Vichnevetsky,  D.  Knight, 
and  G.  Richter,  Eds.  Advances  in  Computer  Methods  for  Partial  Differential  Equa¬ 
tions  -  VII,  IMACS,  New  Brunswick,  pp.  799-805,  1992. 

14.  R.  Biswas  and  J.E.  Raherty,  ‘‘Parallel  and  Adaptive  Methods  for  Hyperbolic  Partial 
Differential  Equations,”  in  R.  Vichnevetsky,  D.  Knight,  and  G.  Richter,  Eds., 
Advances  in  Computer  Methods  for  Partial  Differential  Equations  -  VII,  IMACS, 
New  Brunswick,  pp.  67-74,  1992. 

15.  K.D.  Devine  and  J.E.  Raherty,  ‘‘Parallel  High-Order  Finite  Element  Methods  for 
Conservation  Laws,”  in  R.  Vichnevetsky,  D.  Knight,  and  G.  Richter,  Eds.,  Advances 
in  Computer  Methods  for  Partial  Differential  Equations  -  VII,  IMACS,  New 
Brunswick,  pp.  202-208,  1992  (with  K.  D.  Devine). 

16.  Y.J.  Wang  and  J.E.  Raherty,  ‘‘An  Adaptive  hp-Refinement  Finite  Element  Method 
for  Two-Dimensional  PDEs,”  in  R.  Vichnevetsky,  D.  Knight,  and  G.  Richter,  Eds., 
Advances  in  Computer  Methods  for  Partial  Differential  Equations  -  VII,  IMACS, 
New  Brunswick,  pp.  788-793,  1992. 

17.  M.  Benantar,  J.E.  Raherty,  and  M.S.  Krishnamoorthy,  ‘‘Coloring  Procedures  for  Fin¬ 
ite  Element  Computation  on  Shared-Memory  Parallel  Computers,”  in  A.K-  Noor, 
Ed.,  Adaptive,  Multilevel,  and  Hierarchical  Computational  Strategies,  AMD  Vol, 
157,  ASME,  New  York,  pp.  435-490,  1992. 

Manuscripts  In  Press 

18.  P.K.  Moore  and  J.E.  Raherty,  ‘‘High-Order  Adaptive  Finite  Element-Singly  Implicit 
Runge-Kutta  Methods  for  Parabolic  Differential  Equations,”  BIT,  to  appear,  1993. 

19.  S.  Adjerid,  J.E.  Raherty,  and  Y.J.  Wang,  ‘‘Adaptive  Method-of-Lines  Techniques  for 
Parabolic  Systems,”  Proc.  First  International  Symposium  on  Method  of  Lines.  Sur¬ 
faces  and  Dimensional  Reduction  in  Computational  Mathematics  and  Mechanics, 
1993,  to  appear. 

20.  A.C.M.  Moraes,  J.E.  Raherty,  and  H.  T.  Nagamatsu,  ‘‘Compressible  Laminar  Boun¬ 
dary  Layers  for  Perfect  and  Real  Gases  in  Equilibrium  at  Mach  Numbers  to  30,” 
AIAA  Paper  92-0757,  Thirtieth  Aerospace  Sciences  Meeting  and  Exhibit.  Reno,  Jan. 
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6-9,  1992.  Also,  submitted  for  publication. 

21.  M.  Benantar,  J.E.  Flaherty,  and  M.S.  Krishnamoorihy,  "Triangle  Graphs,"  submitted 
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